#include "fvCFD.H"  
#include "argList.H"  
#include "volFields.H"  
#include "Time.H"  
#include "OSspecific.H"  

int main(int argc, char *argv[])  
{  
    argList::validArgs.append("fieldName");  
    // argList::validOptions.insert("overwrite");  
    #include "setRootCase.H"  
    #include "createTime.H"  
    #include "createMesh.H"  

    word fieldName("T");  
    if(args.args().size() > 1)  
        fieldName = args[1];  

    Info << "Working on field: " << fieldName << endl;  

    volScalarField T  
    (  
        IOobject  
        (  
            fieldName,                // "T"  
            runTime.timeName(),       // e.g. "0", "100"  
            mesh,  
            IOobject::READ_IF_PRESENT,  
            IOobject::AUTO_WRITE  
        ),  
        mesh,  
        dimensionedScalar("zero", dimless, 0.0)  
    ); 

    // 获取整体边界  
    const vectorField& cellCentres = mesh.C();  
    const label nCells = mesh.nCells();  
    
    scalar minX=VGREAT, maxX=-VGREAT, minY=VGREAT, maxY=-VGREAT, minZ=VGREAT, maxZ=-VGREAT;  
    forAll(cellCentres, i)  
    {  
        minX = min(minX, cellCentres[i].x());  
        maxX = max(maxX, cellCentres[i].x());  
        minY = min(minY, cellCentres[i].y());  
        maxY = max(maxY, cellCentres[i].y());  
        minZ = min(minZ, cellCentres[i].z());  
        maxZ = max(maxZ, cellCentres[i].z());  
    }  

    // 四等分的x/y区间  
    scalar midX = (maxX+minX)/2.0, midY = (maxY+minY)/2.0;  
    // 为“字母T”准备粗体行/列宽度  
    scalar xWidth = (maxX-minX)/40.0; // T的横粗，1/40框宽  
    scalar yWidth = (maxY-minY)/40.0; // T的竖粗，1/40框高  
    scalar T_len = (maxY-minY)/8.0; // T的长度，高占1/8  

    // 四个区块中心  
    List<vector> blockCenters(4);  
    blockCenters[0] = vector((minX+midX)/2, (minY+midY)/2, (minZ+maxZ)/2); // 左下  
    blockCenters[1] = vector((minX+midX)/2, (midY+maxY)/2, (minZ+maxZ)/2); // 左上  
    blockCenters[2] = vector((midX+maxX)/2, (minY+midY)/2, (minZ+maxZ)/2); // 右下  
    blockCenters[3] = vector((midX+maxX)/2, (midY+maxY)/2, (minZ+maxZ)/2); // 右上  

    forAll(blockCenters, b)  
    {  
        const scalar cx = blockCenters[b].x();  
        const scalar cy = blockCenters[b].y();  
        const scalar cz = blockCenters[b].z();  

        forAll(cellCentres, i)  
        {  
            const vector& cc = cellCentres[i];  

            // 分块过滤：b=0 左下, b=1 左上, b=2 右下, b=3 右上  
            // bool inX = (b % 2 == 0) ? (cc.x() <  midX) : (cc.x() >= midX);  
            // bool inY = (b < 2)      ? (cc.y() <  midY) : (cc.y() >= midY);  
            // if (!(inX && inY)) continue;  

            // 计算相对中心偏移  
            scalar dx = cc.x() - cx;  
            scalar dy = cc.y() - cy;  

            // 横杆：以块中心为水平居中，厚度为 yWidth，高度方向容差是 yWidth，长度方向半长为 T_len/2  
            bool onTopBar =  
                (Foam::mag(dy) < yWidth) && (Foam::mag(dx) < T_len/2.0);  

            // 竖杆：位于横杆下方（dy < 0），竖向长度 T_len，宽度 xWidth  
            bool onStem =  
                (Foam::mag(dx) < xWidth) && (dy < 0) && (Foam::mag(dy) < T_len);  

            if (onTopBar || onStem)  
            {  
                T[i] = 999.0;  
            }  
        }  
    }  

    // 输出结果  
    T.write();  

    Info << "Field " << fieldName << " with T-shape patch saved." << endl;  

    return 0;  
}